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It is known that quantum computers, if available, would allow an exponential decrease in the 
computational cost of quantum simulations. We extend this result to show that the computation 
of molecular properties (energy derivatives) could also be sped up using quantum computers. We 
provide a quantum algorithm for the numerical evaluation of molecular properties, whose time cost 
is a constant multiple of the time needed to compute the molecular energy, regardless of the size of 
the system. Molecular properties computed with the proposed approach could also be used for the 
optimization of molecular geometries or other properties. For that purpose, we discuss the benefits 
of quantum techniques for Newton's method and Householder methods. Finally, global minima for 
the proposed optimizations can be found using the quantum basin hopper algorithm, which offers 
an additional quadratic reduction in cost over classical multi-start techniques. 



Applying ah initio methods of quantum chemistry 
to particular problems often requires computing deriva- 
tives of the molecular energy. For instance, obtaining a 
molecule's electric properties relies on the ability to com- 
pute derivatives with respect to external electromagnetic 
fields. Likewise, computing the gradient of the molecular 
energy with respect to the nuclear coordinates is the most 
commonly used method for the proper characterization of 
potential energy surfaces and for optimizing the geome- 
try of all but the smallest molecules. The computation of 
these kinds of derivatives, known as molecular properties, 
is nowadays a routine matter when it comes to low-order 
derivatives or small systems (or both). This is largely 
due to advances in analytical gradient techniques, which 
allow for explicit property evaluation without resorting 
to numerical differentiation P dH H 13 16] |7] . 

Nevertheless, the computation of higher-order deriva- 
tives is often prohibitively expensive, even though such 
derivatives are often needed. For example, third- and 
fourth-order anharmonic constants are sometimes re- 
quired to accurately compute a vibrational absorption 
spectrum [3] or efficiently determine the location of tran- 
sition states on complex potential energy surfaces 0. 
Other properties of interest, such as hyperpolarizabilities, 
Raman intensities, or vibrational circular dichroism, are 
also cubic or quartic derivatives. In this report, we show 
that quantum computers, once available, will be able to 
bypass some of the high cost of computing these proper- 
ties. In particular, we show that any molecular property 
can be evaluated on a quantum computer using resources 
that, up to a small constant, are equal to those required 
to compute the molecular energy once. We have previ- 
ously characterized the advantage of quantum computers 
at both computing molecular energies jS] |9] and simulat- 
ing chemical reaction dynamics [lO] , and the present work 
extends our program to molecular properties. 

This paper begins with a brief overview of classical 
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techniques for the evaluation of molecular properties, 
both numerical and analytical. We then introduce the 
quantum algorithm for molecular properties, and discuss 
its advantages and disadvantages with respect to classi- 
cal techniques. We conclude with geometry optimization 
as a particular example, and we show that it can benefit 
from an additional quadratic speed-up through Grover's 
search 1111. 



I. THE CLASSICAL METHODS 

Given an external perturbation /i, the total molecular 
electronic energy can be expanded in a Taylor series 



(1) 



where the coefficients E'"^ are called the molecular prop- 
erties and describe the response of the system to the 
applied perturbation |7]. We consider time-independent 
properties, which can be obtained by differentiating the 
energy at = 0, 
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Many examples of useful derivatives can be given. For in- 
stance, the derivatives with respect to the electric field F 
are the permanent electric dipole, the static polarizabil- 
ity, and the static hyperpolarizabilities of various orders: 
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where the subscript denotes differentiation at F = 0. 
The derivatives with respect to nuclear coordinates R 
include the forces on the nuclei and the force constants, 
while mixed derivatives can provide information such as 
Raman intensities [2]. 

On a classical computer, an energy derivative can be 
evaluated either numerically or analytically, and we dis- 
cuss each approach in turn. 
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Figure 1: Obtaining a numerical gradient of a function de- 
fined on a d-dimensional space classically requires sampling 
the function d + 1 times, once at the origin and once at a 
distance h along each of the axes. Shown above are the sam- 
ple points for the cases d — 1 through d = 3. The quantum 
gradient algorithm can evaluate many sample points in super- 
position, producing the same calculated gradient using only 
one call to the function. 



Numerical derivative techniques rely on computing the 
value of the energy at several discrete points, and then 
using those values to estimate the true derivative. The 
simplest technique is finite difference, which for the first 
derivative in one dimension is the familiar formula, 
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In d dimensions, computing the gradient requires at least 
d+1 evaluations of the energy, once at the origin and once 
at a distance h along each axis (Fig. [ij. Similarly, eval- 
uating higher-order derivatives requires the knowledge of 
the energy on a particular grid, with at least d" + l points 
for the n^^ derivative. 

While numerical gradient techniques usually require 
minimal effort to implement, they are occasionally sus- 
ceptible to numerical instability, due to the ill-posedness 
of numerical differentiation in general [12] . This is partic- 
ularly problematic when using finite-precision arithmetic, 
where various rounding errors can accumulate and be 
amplified upon division by the small number h. The fact 
that small errors in the evaluated function can lead to 
large errors in the derivative affects ab initio electronic 
structure methods insofar as they usually involve long 
calculations with many potential sources of error, includ- 
ing rounding and quadrature. 

By contrast, analytic derivative techniques are those 
that compute the derivative by direct evaluation of an 
analytic expression. They were introduced in quantum 
chemistry by Pulay [1], and have since largely supplanted 
numerical procedures. They are numerically more stable 
and, more importantly, they are usually faster as well. 

Analytic gradient formulas exist for just about all elec- 
tronic structure techniques and for most kinds of pertur- 
bations. To illustrate the argument and establish the cor- 
rect scaling, we will describe the particularly simple case 
of derivatives of fully variational wavefunctions. We start 
by writing the molecular energy as a function E{iJ,, A(/x)) 
of the external perturbation /x and the wavefunction pa- 
rameters A(/x). These parameters, such as the config- 
uration interaction coefficients, completely describe the 
electronic wavefunction. Although A(/x) is a function of 
fi, for simplicity we will write only A. The energy is 



said to be fully variational with respect to A if, for any 
given fi, A assumes the value A* such that the variational 
condition holds: 
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where * indicates A = A*. In that case we can write 
E{(i)=E{(i,X*). 

For fully variational wavefunctions, the gradient with 
respect to fi is given by 
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where we have used the variational condition and the 
Hellman-Feynman theorem. Since one need not know 
the first-order wavefunction response computing the 
gradient is, to within a small constant factor |6], as hard 
as computing the energy. That is, once |A*) is available, 
calculating the expectation value of the Hamiltonian has 
approximately the same computational cost as calculat- 
ing its derivative. However, computing the second deriva- 
tive matrix does require the knowledge of the first-order 
wavefunction response. In fact, as a direct consequence of 
Wigner's 2rt-|- 1 rule of perturbation theory, one needs to 
know the first n wavefunction responses in order to calcu- 
late the (2n -I- 1)*'^ derivative. Computing the responses 
often becomes the bottleneck, and it is what leads to a 
higher asymptotic cost of higher-order derivatives. While 
the gradient requires about the same resources as the en- 
ergy, the second and third derivatives require resources 
that scale as 0{d) times the cost of computing the en- 
ergy (where d is the number of degrees of freedom, i.e., 
the dimension of /i.) |6]. This scaling comes about be- 
cause 0{d) time is required to compute the matrix 

Likewise, the scaling of the (2n-f 1)**^ derivative is 0{d^), 
because the bottleneck becomes the computation of the 
^th Qj.figj. wavefunction response. In other words, the 
computational cost of finding the n*'' analytical deriva- 
tive is 0((iL"/2J)^ roughly a quadratic speed-up over the 
0{d") numerical methods of the same degree. 

The fact that the scaling of derivative techniques, both 
numerical and analytical, depends on d has meant that 
these techniques are often restricted to small systems [29]. 
This is most acutely true of the molecular Hessian, which 
is often beyond reach, even though the gradient is rou- 
tinely accessible. We now show that if quantum com- 
puters were available, the cost of the higher derivatives 
would no longer be prohibitive. 



II. THE QUANTUM ALGORITHM 

The quantum algorithm for molecular properties is 
based on Jordan's quantum gradient estimation algo- 
rithm |1 3) . Jordan's method can numerically compute 
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Table I: Time resources required by various techniques of com- 
puting molecular properties, in terms of the cost of computing 
the energy. For example, the entry "d + 1" means that com- 
puting the property requires d+l evaluations of the molecular 
energy, while the entries in the "Analytical" column indicate 
comparable computational effort. E is the total electronic 
energy, fj, is the external perturbation, and d is the dimen- 
sion of fi. All the derivatives are evaluated at /x = 0. On 
classical computers, the numerical scalings correspond to the 
simplest finite-difference scheme. Analytical techniques are 
the ones that evaluate the derivative directly (the exponent 
[n/2j comes from Wigner's 2n-|- 1 rule). On a quantum com- 
puter, the quantum gradient estimation algorithm is used. It 
should be noted that on a quantum computer, the number of 
evaluations of E needed for any derivative is independent of 
d, and thus of the size of the system. 



the gradient of any function F, given a black box (ora- 
cle) that computes the value of F for an arbitrary input. 
In particular, the algorithm can evaluate the gradient 
using a single query to F, regardless of the number of 
dimensions d of the domain of F. By contrast, the sim- 
plest classical finite-difference scheme would require d-l- 1 
queries to F (see Fig. [T|. The speed-up is essentially 
achieved by being able to sample along all the d dimen- 
sions in superposition. We apply Jordan's algorithm to 
the computation of molecular properties by specifying a 
way to compute the energy on a quantum computer as 
well as by outlining how to obtain higher derivatives. In 
this section, we describe the algorithm, its application to 
quantum chemistry, and finally argue that a return to 
numerical techniques for molecular properties would be 
justified if quantum computers became feasible. 

We assume that the molecular energy is a smooth, 
bounded function of the perturbation, E : 
[^-min 7 £'max] , where a small h is chosen so that E varies 
sufficiently slowly over the domain. For convenience, we 
express the perturbations in units such that h is unitless 
and such that the bounds are the same along all of the 
axes. We also assume that we have a black box for E, 
which, given a perturbation fi, outputs the energy E{fi). 
The precise nature of the algorithm inside the black box 
is not important, so long as it can be implemented on 
a quantum computer. In particular, any classical tech- 
nique of electronic-structure theory can be converted into 
a quantum algorithm 



|14) . In Sec. Ill we will discuss 
possible realizations of the black box, including the use 
of quantum simulation algorithms, which offer a signifi- 



cant improvement over classical ones. 

We begin by choosing the number n of bits of precision 
that we desire in the final gradient. Jordan's algorithm 
starts in an equal superposition on d registers of n qubits 
each {nd qubits total) |14) : 
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where N = 2", the states \ki) are integers on n qubits 
represented in binary, and |k) is a c?-dimensional vector 
of all the |fc)'s. 

We use this state as an input for the black box for E, 
which will, for every integer-vector |k) in the superposi- 
tion, append a phase proportional to the energy E{fj,) at 
perturbation fi — ft,(k — N/2)/iV (where N is the vector 
{N, N, . . . , N) and serves to center the sampling region on 
the origin). To achieve maximum precision with fewest 
qubits, one needs an estimate m of the largest magnitude 
of any of the first derivatives of E. Then, the energy eval- 
uated by the black box is scaled by a factor 2iTN/hm. 
Because the black box operates on all the terms in the 
superposition at once, a single call results in the state 
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The neglect of terms quadratic in h and higher is a valid 
approximation for sufficiently small h (the error caused 
by higher-order terms is discussed in [13] and is only poly- 
nomial). The final state is separable and equals 
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Applying the inverse quantum Fourier transform jl4] 
to each of the d registers results in the gradient 
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The scaling factor N/ m ensures that ^ 3^ is an integer- 
vector with n bits of precision along each dimension. It 
should be reiterated that a single call to E was made, as 
opposed to the d+l that would be needed in the case of 
numerical differentiation by finite difference. 
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Overall, the gradient estimation algorithm produces 
the transformation 



i$(0) 



N dE 
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We can compute the Hessian (and higher derivatives) by 
iterating this algorithm. If, instead of making a call to 
E{fi), the algorithm sought E{fi — i/) from the black box, 
we would perform, at the cost of this single additional 
subtraction, 



N dE 
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Because we will be using this procedure as a subroutine, 
it is important to remove (or "uncompute") the global 
phase, which would otherwise become a relative phase. 
One additional evaluation of E (at i/) suffices for this 
uncomputation. Overall, this supplies another black box. 



which, given i/, yields 



N dE 
m d/x 



using only two calls to 



the original black box for E. One can use the gradient 
algorithm with this new black box, producing the state 
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which is a two-dimensional array of (P quantum registers 
containing all the elements of the Hessian matrix of E. 
In addition, m^^^ is an estimate for the magnitude of the 
largest second derivative, and the phase is 
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Computing higher derivatives would require additional 
factors of two in the number of required black box calls, 
caused by the need to uncompute a global phase at each 
step (this problem is a common feature when it comes 
to recursing quantum algorithms [15]). Hence, evaluat- 
ing the n*'' derivative requires 2"~^ queries to E, which, 
although exponential in n, is much better than d" -I- 1, 
which is the minimum number of function queries re- 
quired to compute the derivative by classical finite dif- 
ference. We stress that the number of calls to E is inde- 
pendent of d, and thus of the size of the system, for the 
derivative of any order. 

One could remark that the quantum gradient algo- 
rithm is a numerical approach and that therefore, just 
like classical numerical techniques, it would be affected 
by numerical instability. This implies that the quantum 
gradient algorithm cannot be used indiscriminately for 
problems that feature errors that cannot be controllably 



reduced through additional computational effort. For ex- 
ample, finite-precision arithmetic presents the same prob- 
lems to quantum computers as it does to classical ones, 
but the rounding errors can be brought under control by 
using more digits of precision (as on classical computers). 
Quantum chemistry techniques might present numerical 
problems as well, insofar as they contain uncontrolled 
sources of error. However, if the technique for comput- 
ing the energy is numerically exact, that is to say, if the 
error in the energy can be controllably reduced below any 
level, the magnitude of the numerical error in the calcu- 
lated derivative can likewise be made arbitrarily small. 
Fortunately, quantum computers would make it possible 
to efficiently evaluate the numerically exact molecular 
energy, meaning that numerical instability will not be a 
problem. We turn to the topic of molecular energy eval- 
uation next. 



III. THE BLACK BOX FOR THE ENERGY 

The application of Jordan's gradient algorithm to 
chemical problems requires that there be a black box 
that can compute the value of the ground-state molec- 
ular energy at any value of the perturbation /x in the 
neighborhood of /x = 0. Furthermore, to avoid numeri- 
cal artifacts, this black box should be numerically exact, 
allowing the error in the energy to be controllably re- 
duced through additional computational work. 

The problem of exact classical electronic structure 
methods is that they generally have a computational cost 
that scales exponentially with the size of the system. Al- 
though these classical algorithms could also be used as 
subroutines in the quantum gradient algorithm, there are 
quantum electronic-structure algorithms that could avoid 
the exponential scaling in many cases. 

In particular, we have recently described a quantum 
full CI algorithm [8] for computing the molecular ground 
state energy to a given precision in O (M^) time [16], 
where M is the number of basis functions. This algorithm 
could be easily recast as a subroutine that would function 
as the black box for the energy. Several modifications 
would have to be made, including the direct computa- 
tion of the overlap integrals on the quantum computer, 
rendering it possible to introduce the perturbation /j, into 
the calculation. Nevertheless, a quantum computer run- 
ning the quantum FCI algorithm could be used to obtain 
a molecular property of a system with basis size M in 
O (M^) time, a dramatic improvement over the possibil- 
ities of classical computers. 

A more recent development is the real-space chemical 
dynamics simulation algorithm |10| 117). based on Zalka's 
earlier work [18]. It is known that simulating, to a given 
precision, the exact dynamics of a system of P parti- 
cles interacting under a pairwise interaction requires at 
most 0{P^) time and 0{P) space, in contrast to the 
classical exponential cost. If an eigenstate of the sys- 
tem Hamiltonian were prepared as the initial state |19j . 
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the dynamics would only apply a phase to the wavefunc- 
tion. This phase could be read out by the phase estima- 
tion algorithm |14l |20] [2T| , forming the required energy 
black box. Although the large prefactor of this algorithm 
would make it slower for small molecules than the equiv- 
alent quantum FCI calculation, it benefits from a supe- 
rior asymptotic scaling as well as from the fact that only 
minimal modifications would need to be made to insert 
the perturbation fi into the calculation. For example, 
simulations with different nuclear coordinates proceed in 
exactly the same way, while an electromagnetic field re- 
quires only a small modification of the simulated Hamil- 
tonian 

It should be remarked that current quantum algo- 
rithms for energy estimation, such as the ones mentioned 
above, rely on quantum phase estimation, which has been 
criticized as inefficient |2l] because its cost grows expo- 
nentially with the number of bits of precision sought. 
This could be significant for gradient estimation, which 
may require precise energy evaluations to avoid numerical 
errors. To estimate the cost, we note that if the gradient 
is desired to n bits of precision (as in Eq. Ill, the black 
box should evaluate the energy to 
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bits of precision [13] , where cos^ 9 is the desired success 
probability of the algorithm. For example, with 6 = tt/8, 
the algorithm succeeds 85% of the time and requires four 
more digits of precision in the energy than is desired in 
the gradient. The four additional digits present only 
a constant overhead, meaning that the computation of 
any molecular property at any precision is, up to a con- 
stant factor, as hard as computing the energy of the same 
molecule at the same precision. 

Finally, a limitation of current quantum simulation al- 
gorithms is that they are generally spin-free and non- 
relativistic, which limits the ability to compute deriva- 
tives such as indirect spin-spin coupling. 



IV. NEWTON'S METHOD AND GEOMETRY 
OPTIMIZATION 

Perhaps the single most common use of molecular 
derivatives is molecular geometry optimization. We can 
therefore use it to illustrate some of the advantages of 
a quantum algorithm over a classical one, including a 
quantum version of Newton's method, which offers an 
additional quadratic speedup over its classical counter- 
part. 

A simple way for finding the locally optimal geometry 
is to perform the standard Newton iterations. 



dR2 



dR2 
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until convergence is reached. Here, R„ are the nuclear 
coordinates at the n*'^ iteration, and and 

are, respectively, the gradient and Hessian of E with re- 
spect to nuclear displacement (the "molecular gradient" 
and the "molecular Hessian"). If a quantum computer 
were used to compute the derivatives, one would require 
exactly 3 calls to a black box for E per iteration: one 
for the gradient and two for the Hessian. A classical ap- 
proach, on the other hand, would be much slower, requir- 
ing at least (P + 1 function calls for finite difference, and 
approximately 0{d) effort in the analytical case j30]. For 
large molecules with large d, this savings could prove sig- 
nificant, even if each energy evaluation takes much longer 
on a quantum machine than on a classical computer. 

There are many classical tricks available for speed- 
ing up the convergence of Newton's method if the initial 
guess is not close to a local minimum, in which case the 
usual Newton step might be inappropriately large. Tech- 
niques such as trust regions and level shifts [23] are still 
available to quantum computers, or they can be imple- 
mented as classical post-processing. 

In addition, we remark that Newton's method is the 
first in the class of Householder methods, which offer a 
rate of convergence of £ + 1, provided that derivatives 
up to order £ + 1 exist and can be calculated. A quan- 
tum computer could be used to accelerate Householder 
methods of any degree, requiring 2™"^ — 2^+^ — 1 

calls to the black box for order-£ Householder optimiza- 
tion method. Although exponential in £, this expression 
is independent of system dimension d. 

Of course, Newton's method is only useful for local 
minimization, and we are often interested in global op- 
timization. Here, we can use a quantum version of the 
multistart technique, called the quantum basin hopper 
|24l l25l l26l [27] . A number of points is selected at ran- 
dom, and each is followed, using a local search, to its local 
basin (if a quantum version of Newton's method is used 
for the local search, such as the one we propose above, we 
can get the usual quadratic convergence). Then, the min- 
ima of all the basins are compared and the least one cho- 
sen as the global minimum. Quantum computers could 
add a quadratic speed-up to such a multistart technique, 
since the resulting local minima form an unstructured 
database that can be searched using Grover's algorithm 
jlll [T3] with a quadratic speed-up. As Diirr and H0yer 
pointed out |28], a Grover search can find the minimum 
of an unstructured database with 0{^/K log K) calls to 
the database (where K is the number of database en- 
tries, i.e. multistart points), as opposed to the classically 
required 0{K log K) queries. 



V. CONCLUSION 

We have shown that Jordan's quantum gradient es- 
timation algorithm can be applied to the estimation of 
time-independent, non-relativistic molecular properties. 
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Doing so requires a quantum electronic-structre black 
box, for which known quantum simulations algorithms 
are well suited. The quantum algorithm offers a speed- 
up from the classical cost of 0(dL"/2J) ^^-^ analj'tical 
derivatives to the quantum query complexity of 2"'~^. 
That is, the number of energy calculations required on 
the quantum computer is independent of d, and thus of 
the system size, which could offer a significant advantage 
for the computation of properties of large systems. In 
particular, it would make the molecular Hessian of any 
molecule only twice as expensive as its molecular gra- 
dient, enabling a fast, local geometry optimization us- 



ing Newton's method. Finally, global optimization could 
combine the local Newton's method with Grover search 
to offer an additional quadratic speed-up over the classi- 
cal multi-start technique. 
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